/*
 * do_rdf.c
 *
 *  Created on: Apr 14, 2010
 *      Author: marchi
 */

#include "stdio.h"
#include "stdlib.h"
#include "statutil.h"
#include "typedefs.h"
#include "smalloc.h"
#include "vec.h"
#include "copyrite.h"
#include "statutil.h"
#include "math.h"
#include "rdf_grid.h"
#include "xvgr.h"


int do_rdf_grid(real*** Grid, int nx, int ny, int nz, matrix co, matrix oc, rvec com, hist_t *hgrid){
//	FILE *fp;
//	fp=fopen("Mydata.dat","w+");
	int i,j,k,h;
	rvec xa,com1,x,xc;
	real dist;
	com1[XX]=oc[XX][XX]*com[XX]+oc[XX][YY]*com[YY]+oc[XX][ZZ]*com[ZZ];
	com1[YY]=                   oc[YY][YY]*com[YY]+oc[YY][ZZ]*com[ZZ];
	com1[ZZ]=									   oc[ZZ][ZZ]*com[ZZ];
	for(i=0;i<nx;i++){
		xa[XX]=(real)i/(real)nx;
		for(j=0;j<ny;j++){
			xa[YY]=(real)j/(real)ny;
			for(k=0;k<nz;k++){
				xa[ZZ]=(real)k/(real)nz;
				x[XX]=xa[XX]-com1[XX];
				x[YY]=xa[YY]-com1[YY];
				x[ZZ]=xa[ZZ]-com1[ZZ];
				x[XX]=x[XX]-rint(x[XX]);
				x[YY]=x[YY]-rint(x[YY]);
				x[ZZ]=x[ZZ]-rint(x[ZZ]);

				xc[XX]=co[XX][XX]*x[XX]+co[XX][YY]*x[YY]+co[XX][ZZ]*x[ZZ];
				xc[YY]=                 co[YY][YY]*x[YY]+co[YY][ZZ]*x[ZZ];
				xc[ZZ]=                                  co[ZZ][ZZ]*x[ZZ];
				dist=sqrt(sqr(xc[XX])+sqr(xc[YY])+sqr(xc[ZZ]));
				if(dist<hgrid->rmax){
					h=(int) (dist/hgrid->binsize);
					if(h > hgrid->nmax) gmx_fatal(FARGS,"stop in the middle");
					hgrid->hist[h]+=Grid[i][j][k];
					hgrid->nhist[h]++;
				}
			}
		}
	}
//	gmx_fatal(FARGS,"There are no molecule descriptions. I need a tpr file for this pbc option.");
	return 0;
}

int setup_rdf_grid(hist_t *hgrid){
	hgrid->binsize=0.005;
	hgrid->rmax=2.0;
	hgrid->nmax=(int) ceil(hgrid->rmax/hgrid->binsize);
	snew(hgrid->hist,hgrid->nmax);
	snew(hgrid->nhist,hgrid->nmax);
	return 0;
}

int print_rdf_grid(FILE *fp,const char *distfile, hist_t hgrid, const output_env_t oenv){
	int i;
	real dx,x,fx;
	dx=hgrid.binsize;
	for(i=0;i<hgrid.nmax;i++){
		x=dx*((real) i);
		if(hgrid.nhist[i]) fx=hgrid.hist[i]/((real)hgrid.nhist[i]);
		else fx=0.0;
		if(fx != 0.0) fprintf(fp," %g  %g\n",x,(double) fx);
	}
	return 0;
}

int do_rdf_dip_grid(real*** Grid, int nx, int ny, int nz, matrix co, matrix oc, rvec com, rvec dip, hist_t *hgrid){
//	FILE *fp;
//	fp=fopen("Mydata.dat","w+");
	rvec dip_n;

	int i,j,k,h;
	rvec xa,com1,x,xc;
	real dist, dip_mag;
	unitv(dip,dip_n);
	dip_mag=norm(dip);
	com1[XX]=oc[XX][XX]*com[XX]+oc[XX][YY]*com[YY]+oc[XX][ZZ]*com[ZZ];
	com1[YY]=                   oc[YY][YY]*com[YY]+oc[YY][ZZ]*com[ZZ];
	com1[ZZ]=									   oc[ZZ][ZZ]*com[ZZ];
	for(i=0;i<nx;i++){
		xa[XX]=(real)i/(real)nx;
		for(j=0;j<ny;j++){
			xa[YY]=(real)j/(real)ny;
			for(k=0;k<nz;k++){
				xa[ZZ]=(real)k/(real)nz;
				x[XX]=xa[XX]-com1[XX];
				x[YY]=xa[YY]-com1[YY];
				x[ZZ]=xa[ZZ]-com1[ZZ];
				x[XX]=x[XX]-rint(x[XX]);
				x[YY]=x[YY]-rint(x[YY]);
				x[ZZ]=x[ZZ]-rint(x[ZZ]);
				xc[XX]=co[XX][XX]*x[XX]+co[XX][YY]*x[YY]+co[XX][ZZ]*x[ZZ];
				xc[YY]=                 co[YY][YY]*x[YY]+co[YY][ZZ]*x[ZZ];
				xc[ZZ]=                                  co[ZZ][ZZ]*x[ZZ];
				rvec xd;
				unitv(xc,xd);
				real theta=gmx_angle(xd,dip_n);
				dist=sqrt(sqr(xc[XX])+sqr(xc[YY])+sqr(xc[ZZ]));
				if(dist<hgrid->rmax && abs(RAD2DEG*theta-90.0) > 30.0){
					h=(int) (dist/hgrid->binsize);
					if(h > hgrid->nmax) gmx_fatal(FARGS,"stop in the middle");
					hgrid->hist[h]+=Grid[i][j][k]/cos(theta)/dip_mag;
					hgrid->nhist[h]++;
				}
			}
		}
	}

//	gmx_fatal(FARGS,"There are no molecule descriptions. I need a tpr file for this pbc option.");
	return 0;
}

